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Time-dependent modeling of pulsar wind nebulae: 
Study on the impact of the diffusion-loss approximations 
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ABSTRACT 

In this work, we present a leptonic, time-dependent model of pulsar wind nebulae 
(PWNe). The model seeks a solution for the lepton distribution function consider- 
ing the full time-energy dependent diffusion-loss equation. The time-dependent lep- 
ton population is balanced by injection, energy losses, and escape. We include syn- 
chrotron, inverse Compton (IC, with the cosmic- microwave background as well as 
with IR/optical photon fields), self-synchrotron Compton (SSC), and bremsstrahlung 
processes, all devoid of any radiative approximations. With this model in place we 
focus on the Crab nebula as an example and present its time dependent evolution. 
Afterwards, we analyze the impact of different approximations made at the level of 
the diffusion-loss equation, as can be found in the literature. Whereas previous mod- 
els ignored the escape term, e.g., with the diffusion- loss equation becoming advective, 
others approximated the losses as catastrophic, so that the equation has only time 
derivatives. Additional approximations are also described and computed. We show 
which is the impact of these approaches in the determination of the PWN evolution. 
In particular, we find the time-dependent deviation of the multi-wavelength spectrum 
and the best-fit parameters obtained with the complete and the approximate models. 

Key words: pulsars: general, radiation mechanisms: non-thermal 



1 INTRODUCTION 

In addition to their electromagnetic emission, pulsars dis- 
sipate their rotational energy via relativistic winds of par- 
ticles. Because the relativistic bulk velocity of the wind is 
supersonic with respect to the ambient medium, such a wind 
produces a termination shock. In turn, the wind particles, 
moving trough the magnetic field and the ambient photons, 
produce radiation that we observe as Pulsar Wind Nebulae 
(PWNe). As the pulsars themselves, the PWN emits at all 
wavelengths from radio to TeV. 

PWNe usually have two main X-ray morphologies, de- 
pending to the velocity of the pulsar proper motion. For slow 
pulsars, images taken with the Chandra X-ray Observatory 
(see e.g., Kargaltsev & Pavlov, 2008) show a toroidal shape 
around the pulsar equator, with two possible jets starting 
from the pulsars poles. Instead, pulsars moving with high 
velocity in the interstellar medium produce PWNe with the 
characteristic bullet-like or bow-shock morphology, with the 
tail developed along the pulsar motion. Thus, the study of 
PWNe can lead to knowledge of pulsar winds, the properties 
of the ambient medium, and the wind-medium interaction. 

PWNe constitute the largest class of identified Galac- 
tic very-high energy (VHE) gamma-ray sources, with the 
number of TeV detected objects increasing from 1 to ~30 



in the last 6 years. This statistics shines in comparison with 
the ~30, 10, or 40 PWNe known in radio, optical/IR, or X- 
rays, respectively, detected in decades of observations. The 
majority of PWNe at VHE gamma-rays have a very large 
size, depending on their evolutionary stage and proximity 
(see e.g., Vela X or HESS J1825-137). The Cherenkov Tele- 
scope Array (CTA, Actis et al. 2012), with its large FoV, 
will be able to image the whole plerionic emission, including 
the halo produced by cooled electrons. The larger energy 
range of CTA compared with all other gamma-ray instru- 
ments will be crucial to understand cooling effects, includ- 
ing the resolution of internal structures and the properties 
which could be therein active, such as their likely enhanced 
magnetic field, and the possibility of disentangling between 
synchrotron and adiabatic losses. For the first time, we will 
be able to study in detail the interaction between the host 
remnant and its PWN, and how much of the gamma-ray 
emission originates in each component should feedback the 
theoretical understanding of the evolutionary tracks. Ob- 
servations with CTA will hopefully produce a homogeneous 
sampling of the Galactic PWNe, since its sensitivity will 
permit the detection of PWNe up to 50 kpc. Out of the 
complete PWNe population in the Galaxy, assuming 40 kyr 
for the estimated lifetime of TeV-emitting electrons for a 
magnetic field of 3 micro Gauss, CTA would detect 300-600 



2 Martin, Torres, Rea 



objects, of which between 15-25% will be fully resolved, de- 
pending on proximity, age, and flux (see de Ona et al. 2012 
for detailed studies on CTA expectations). 

In studying PWNe, there are two distinct theoret- 
ical approaches. On the one hand, detailed magneto- 
hydrodynamic (MHD) simulations have succeeded in ex- 
plaining the morphology of PWNe (see Rea & Torres 2011 
for reviews on several aspects of these issues). On the other 
hand, spherically symmetric ID PWNe spectral models, 
with no energy-dependent morphological output, have been 
constructed since decades (although there are only a handful 
of such codes, and none is public to our knowledge), see, e.g., 
Aharonian et al. 1997, Atoyan & Aharonian 1998, Bednarek 
& Bartosik 2003, ibid. 2005; Biisching et al. 2008; Fang & 
Zhang 2010; Zhang et al. 2008; Tanaka & Takahara 2010, 
2011; Li et al. 2010; and other references quoted below. 

Some of the most elaborate models solve the time- 
energy dependent diffusion-loss equation, where the time- 
dependent lepton population is balanced by injection, energy 
losses, and escape, with various degrees of approximations. 
Such models of PWNe usually calibrate with the Crab neb- 
ula measurements today, being it the best studied PWN 
at all wavelengths. Whereas essentially all models prop- 
erly fit the Crab nebula data at the current age, the time- 
evolution out of the normalization age of the PWN can show 
significant deviations. In this work we present a leptonic, 
time-dependent model of pulsar wind nebulae (PWNe). It 
seeks a solution for the electron distribution function con- 
sidering the full time-energy dependent diffusion-loss equa- 
tion. The time-dependent lepton population is balanced 
in our model by injection, energy losses, and escape. We 
include synchrotron, Inverse Compton (with the cosmic- 
microwave background as well as IR/optical photon fields), 
self-synchrotron Compton, and bremsstrahlung processes, 
all devoid of any radiative approximations. With this model 
in place we focus on the Crab nebula as an example and 
achieve a fit for its persistent emission. Afterwards, we an- 
alyze the impact of different approximations made at the 
level of the diffusion-loss equation, as can be found in the 
literature. We note that conclusions on specific sources and 
population analysis using approximate tools is severely af- 
fected. In §2 we describe our model. In §3, we show the 
results obtained for the Crab nebula and the luminosity ra- 
tios computed by our code and compared with observations. 
In §4, we comment on different approximations found in the 
literature and we show their impact in the Crab nebula evo- 
lution. Finally, in §5, we write our concluding remarks. 



2 MODELING THE EMISSION OF PWNE 
2.1 The diffusion equation 

The diffusion equation adopted in this work reads (e.g. 
Ginzburg & Syrovatsky 1964) 



dt 



-|- W7 , t) iV( 7) *)]-fe| + Q( 7l i), (1) 

where the left-hand side is the variation of the lepton dis- 
tribution in time and the first term in the right-hand side 
accounts for the continuous change in energy of the particles 
due to energy losses. The function 7(7, t) is the summation 
of the energy losses due to all processes considered. Q(7, t) 



represents the injection of particles per unit energy and unit 
volume in a certain time. r( 7 , t) represents the escape time, 
after which the particles are effectively removed from the 
phase space. We assume that particles escape (second term 
in the right hand side of Eq. 1) via Bohm diffusion (e.g., as 
in Zhang et al. 2008, or Li et al. 2010). Eq. ([TJ is solved using 
the implicit forward difference technique on the derivatives 
in time and energy. 

2.2 The injection of particles 

Our numerical implementation allows for any form of parti- 
cle injection. We assume a broken power-law 



Q(7,t) =Qo(t) 



7b / 



for 7 < 7b, 

(2) 

for 7 > 76, 



where 7;, is the break energy. The parameters ai and Q2 
are the spectral indices. The normalization constant Qo(t) 
is determined using the injection luminosity L(t), 



Lit) = L 1 + — 

To 



(3) 



and where Lo is the initial luminosity, to is the initial spin- 
down timescale of the pulsar and n is the breaking index. 
These parameters may be observationally determined (see, 
e.g., Lewin & van der Klis 2006, Gaensler & Slane 2006). In 
the case of the spin-down luminosity, we have 



L{t )=4n 2 I^, 



(4) 



where P and P are the period and its first derivative and 
/ is the pulsar moment of inertia, which we assume I ~ 
10 45 g cm 2 . The initial spin-down timescale of the pulsar is 
(Gaensler & Slane 2006) 

Po 2r c 



To 



(n - 1)P n 



1 



tag 



(5) 



(6) 



where Po and Po are the initial period and its first derivative 
and t c is the characteristic age of the pulsar, 

_ P_ 

Tc ~ 2P' 

The breaking index can also be computed from observational 
data when the second derivative of the period, P, is known. 
Assuming that the angular frequency f2 = 2n/P of the pul- 
sar evolves in time as Cl = fcO" where n is again the breaking 
index and k is a constant that depends on the magnetic mo- 
ment of the pulsar, we find n — QQ./Q, 2 ~ PP/P 2 . If the 
system is a dipole spin-down rotator, the breaking index is 
exactly 3 and the constant k has the value k = 2fj 2 ± /3Ic 3 , 
where fi± is the component of the magnetic dipole moment 
orthogonal to the rotation axis. 

The normalization of the injection function is given by 



poo 

(l- ?7 )L(t)= / 7 mc 2 Q( 7 ,t)d7, 
Jo 



(7) 



where rj is the magnetic energy fraction. It is defined as 
77 = Ls(t)/L(f), where Ls(t) is the magnetic power; thus, 
7] is its ratio with the spin-down power. This definition, see 
e.g., Tanaka & Takahara (2010), divides the energy injection 
from the pulsar into magnetic field energy and relativistic 
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particle energy and is different from the one used for the 
magnetization parameter a(t) = Lpit) / L v {t), and where 
L p (t) is the relativistic particle's fraction of the spin-down 
power. To ensure particle confinement, we impose that the 
Larmor radius of the particles has to be smaller than the 
termination shock radius, what leads to 



''/max (£) 



L(t) 



(8) 



where e is the electron charge and e is the fractional size of 
the radius of the shock, which has to be smaller than 1. k is 
the magnetic compression ratio. For strong shocks (a <C 1), 
k ~ 3 (Venter & de Jager 2006, Holler et al. 2012). 

Simulations with relativistic shocks in unmagnetized 
plasmas predict a particle spectrum downstream that has 
two components: a relativistic Maxwellian and a high-energy 
tail fitted by a power-law with an energy index of —2.4 ±0.1 
(Spitkovsky 2008). Some papers extrapolate these results to 
the case of PWNe, see, for instance, Grondin et al. 2011 and 
Van Etten & Romani 2011. However, Spitkovsky's simula- 
tions of relativistic shocks in unmagnetized plasmas acceler- 
ate particles until 7 « 1000 only. Here, for simplicity and to 
facilitate comparison with other models (see e.g., Zhang et 
al. 2008, Gelfand et al. 2009 or Tanaka & Takahara 2010), we 
consider a pure broken power-law injection shown in Eq. ((2j. 

2.3 Energy losses and photon spectrum 

We consider synchrotron, (Klein-Nishina) inverse Compton, 
bremsstrahlung, and adiabatic losses (and their time depen- 
dence). Details of their implementation are given in the Ap- 
pendix. 



3 RESULTS 

3.1 The Crab nebula 

The distance to the Crab Nebula is 2 kpc (Manchester et al. 
2005). The period and its derivative are obtained from Tay- 
lor et al. (1993). Assuming that the moment of inertia of the 
Crab pulsar is I = 10 45 g cm 2 , and using Eq. @, we obtain 
the spin-down luminosity power today. The expansion of the 
PWN is considered using the free expansion approximation 
given by van der Swaluw (2001). We consider a characteris- 
tic energy for the SN explosion of 10 erg and an ejected 
mass of 9.5M@ (Bucciantini et al. 2011). All this parame- 
ters used for the Crab Nebula are summarized in Table [T] 
including those determined by the code. 

At the date of the pulsar period ephemeridis (year 
1994), the age of the pulsar was 940 yr. This is consistent 
with Eq. Q and helps to minimize the bias produced by the 
non-simultaneity of the multi-wavelength data points used, 
obtained from ~1970 (radio) to 2008 (VHE). We checked 
that changing the ephemeris to the latest one (e.g., the one 
used by the Fermi-LAT CollaboratiorQ) introduce no visible 
change in the results. 

In order to compute the IC energy losses and spectrum, 
we consider three components: the cosmic microwave back- 
ground (CMB), the galactic far infrared background (FIR) 

1 http : //f ermi . gsfc.nasa.gov/ssc/data/access/lat/ephems/ 
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Figure 2. Cooling times at t age = 940 yr. At low energies, the 
adiabatic losses are dominant because their cooling time is of the 
same order of the pulsar age. At high energies, synchrotron losses 
become the most important. 



and the near infrared and optical photon field due to the 
stars (NIR), n{v) = uomb{v) + n F m(^) + n NIR (v). Each 
one of the latter two is considered as a diluted blackbody 
(Schlickeiser 2002). The temperature of the FIR (NIR) is 
considered as 70 (5000) K. The CMB is a blackbody of tem- 
perature 2.73 K. 

We find different ways to evolve the magnetic field in the 
literature (see Rees & Gunn 1974, Kennel & Coroniti 1984, 
Venter et al. 2006, de Jager et al. 2009, Tanaka & Takahara 
2010,2011). We assume magnetic energy conservation as in 
Tanaka & Takahara 2010, 



r)L(t')dt' = —RpwN(t) 



47T„ 3 t ^B 2 {i) 



(9) 



thus, using Eq. (|3]) and solving for the field we obtain 



B(t) 



3(w — 1)t?Lqto 



1 + 



(10) 



The magnetic energy is conserved in the sense that the mag- 
netic fraction (the fraction of the spin down power that goes 
into the magnetic field) is constant. This parameter is called 
rj in Eq. (10). The magnetic field itself is thus time depen- 
dent, and its behavior is given in Eq. (|10p . This approach 
has a very similar behaviour as it is adopted in other papers 
cited before. 

For the injection we use Eq. ([2}, where Qo(t) is calcu- 
lated using Eq. (|7|). The final luminosity power is given by 
Eq. Q and the initial spin-down power is determined us- 
ing Eq. (J3j) , since we know the luminosity power nowadays 
and the age of the pulsar. For the ISM density in the Crab 
Nebula, we take a fiducial value of 1 cm" 3 . Thus, our free 
parameters in order to fit the spectrum are the magnetic 
fraction 77 and the shock radius fraction e. For the final fit, 
we get 77 = 0.012 and e = 1/3. The magnetic field value we 
get today is 97 /jG, which is close to the 100 fiG calculated 
by the MHD simulations done by Volpi et al. (2008). Ta- 
ble Q] clarifies which parameters comes from observations or 
assumptions, and which parameters are used to fit the data. 
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Table 1. Summary of the physical magnitudes used or obtained for the Crab Nebula fit at the current age. A few parameters are fixed 
based on prior input or hypothesis. 



Magnitude 



Symbol 



Origin or Result 



A 

Age (yr) 


tage 


940 




Period fmsl 






trom laylor ct al. ^lyyoj 


Period derivative (s s "*") 


P(t a ge) 


4.209599 Xl0~ 13 


from Taylor ct al. (1993) 


Spin-down luminosity now (crg/s) 


L(t age ) 


4.53 X 10 38 


Eq. {3j 


Moment of inertia (g cm 2 ) 


I 


10 45 


Eq. {4j 


Breaking index 


n 


2.509 


from Lyne et al. (1988) 


Distance (kpc) 


d 


2 


from Manchester ct al. (2005) 


Ejected mass (Mq) 




9.5 


from Bucciantini ct al. (2011) 


SN explosion energy (erg) 


Eq 


10 51 


from Bucciantini ct al. (2011) 


Minimum energy at injection 


Imin 


1 


assumed 


Maximum energy at injection at t age 


7m ax {tage ) 


7.9 X 10 9 


result of the fit 


Break energy 


7b 


7 x 10 5 


result of the fit 


Low energy index 


"I 


1.5 


result of the fit 


High energy index 


OC2 


2.5 


result of the fit 


Shock radius fraction 


e 


1/3 


result of the fit 


Initial spin-down luminosity (crg/s) 


La 


3.1 x 10 39 


result of the fit 


Initial spin-down age (yr) 


T o 


730 


Eq. {5j 


Magnetic field (^G) 


S(tage) 


97 


result of the fit 


Magnetic fraction 


V 


0.012 


result of the fit 


PWN radius today (pc) 


RPW N^age) 


2.1 


Eq. f^2l 


CMB temperature (K) 


T CMB 


2.73 


fixed 


CMB energy density (eV/cm^) 


V>CMB 


0.25 


fixed 


FIR temperature (K) 


Tfir 


70 


as in Marsden et al. (1984) and subsequent refs. 


FIR energy density (eV/cm^) 


U>FIR 


0.5 


as in Marsden et al. (1984) and subsequent refs. 


NIR temperature (K) 


Tnir 


5000 


as in Aharonian et al. (1997) and subsequent refs. 


NIR energy density (eV/cm^) 


U>NIR 


1 


as in Aharonian et al. (1997) and subsequent refs. 


Hydrogen density (cm ^) 


Iff 


1 


assumed 
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Figure 1. From top to bottom, left to right: Magnetic field, spin-down power, lepton population, and spectral energy distribution of the 
Crab nebula as a function of time. 
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3.2 The Crab fitting 

Figure [T] shows the magnetic field, spin-down power, lep- 
ton population, and spectral energy distribution of the Crab 
nebula as a function of time, resulting from our code after 
normalization to current measurements. The current cool- 
ing times for the different processes considered are shown in 
Figure[5] whereas the current spectrum is shown in Figure[3] 

The SSC flux is the strongest contributor to the high- 
energy spectra, followed by IC with the CMB and the FIR. 
The Bremsstrahlung contribution is not very important, but 
as it is similar to the NIR radiation, we do not neglect it 
in favor of the other contributions. Most of the radiative 
considerations of Tanaka & Takahara (2010) are similarly 
obtained in our model, since they are driven by SSC dom- 
ination. Our resulting value of the magnetic field today is 
lower than that used by Atoyan & Aharonian (1996) in their 
time-independent approach, who in turn adopted it from the 
Kennel & Coroniti (1984) model, followed by an adjustment 
on the relativistic particle density to enable the data fitting. 
This value of magnetic field is unrealistically high for our 
time-dependent spectral model, and a lower value is pre- 
ferred also by MHD simulations. 

Regarding the time evolution presented in Figure [T] it 
is interesting to note how the peak of the electron distribu- 
tion moves from lower Lorentz factors to the energy break 
in the injection. This displacement of the peak is due to the 
high energy losses for energies lower than the break at early 
ages. The maximum energy of the injection is decreasing 
with time and the maximum energy of the electrons popula- 
tion is decreasing also through energy losses, but at a slower 
rate due to the presence of high-energy electrons that were 
previously injected. The slope of the distribution at VHE 
becomes flattened with time also due to evolution in time of 
the dominant cooling process, increasing the power of the IC 
radiation. As the magnetic field falls, the synchrotron radi- 
ation diminishes with respect the IC radiation and at later 
ages (e.g., towards 10 kyr), the IC radiation contains most 
of the emitted flux. This is in agreement with the idea of 
older PWNe being still detectable at high energies but being 
devoid of lower-energy counterparts (de Jager & Djannati- 
Ata'i 2008). This is shown in Figure 0] We see that the flux 
at energies >1 TeV and the gamma-ray flux are equal for 
an age of ~5 kyr. 

The radio and optical evolution of the Crab nebula show 
a decreasing-with-time behavior. Measurements of the radio 
flux decrease were done by Vinyaikin (2007), using data from 
1977 to 2000 at 86, 151.5, 927 and 8000 MHz. The mean 
flux-decrease rate averaged obtained was -0.17±0.02% yr _1 . 
Using data obtained from our code at the same frequencies 
for the same time interval we obtained an averaged rate 
of -0.2% yr -1 . In optical frequencies, the continuum flux 
decrease 0.5±0.2% yr _1 at 5000 A(Smith 2003). In this case, 
we obtained directly from the model a flux decrease of 0.3% 
yr _1 . The evolution of both luminosities as extracted from 
our model is in agreement with observations. 



4 APPROXIMATE MODELS 

Apart from the approximations we focus below, one can also 
find many radiative approximations too in PWNe models: 



L(>1 TeV)/L(2-10 koV) 

L(0.1-l TcV)/L(3-10 kcV)- 
L(>1 TeV)/L(0.1-l TeV) - - 



1000 
Time (yr) 



Figure 4. Luminosity ratios for the Crab nebula evolution: L(>1 
TeV) / L{2 - 10 KcV), L(0.1 - ITeV) / L(2 - 10 KeV) and L(>1 
TeV) / L(0.1 - 1 KeV) 



using a priori guesses for which field is dominant in each 
environment, using mono-chromatic assumptions for syn- 
chrotron and inverse Compton, or using Thompson cross sec- 
tion instead of Klein-Nishina. These assumptions certainly 
simplify the treatment, but at the expense of assuming ap- 
proximations for which their impact is usually not checked. 
We have not adopted any of them here. 

Regarding the diffussion-loss equation, the most usual 
approximation is to neglect the escape term (see, e.g. Tanaka 
& Takahara 2010, 2011), to obtain an advective differential 
equation (ADE). Using just this approximation in our com- 
plete model would lead to very similar values for the mag- 
netic field and magnetic fraction (needed to obtain a good fit 
for today's Crab nebula, when imposing a correct contribu- 
tion of the SSC such that it fits the high energy data today). 
This is because the Bohm timescale is larger than the age of 
the Crab nebula and is not affecting strongly the particles' 
evolution. Another common (and additional) approximation 
is neglecting the treatment of energy losses and instead re- 
place it by the particle's escape time (see, e.g., Zhang et al. 
2008; Qiao & Fang 2009). In this case, Eq. Q has the form 

dNfrt) _ N( 7 ,t) 



+ Q(7,*). 



(ii) 



dt r(j,t) 

wherer(7,f) = t/ItCTj*)! ' s the escape time of the particles. 
In this case, particles are not losing energy, but they are 
rather removed from the distribution after a certain time. 
This makes Eq. (fTT} a partial differential equation in time 
only (TDE). 

Before doing the fits, we fixed the parameters which 
are obtained from observations, as we have done in the com- 
plete model, and included the ADE and TDE cases with the 
complementary approximations done by Tanaka & Takahara 
(2010, hereafter ADE-T) and Zhang et al. (2008, hereafter 
TDE-Z). In ADE-T, the bremsstrahlung energy losses and 
its spectrum, and the FIR and NIR contributions into the 
inverse Compton energy losses and their spectrum, are ig- 
nored. Also, the maximum energy at injection is fixed and 
the expansion of the PWN is modeled in a ballistic approxi- 
mation (Rpwn = vpwNt). All these approximations are not 
done in the full treatment presented above, against which we 
compare. In the TDE-Z, only the synchrotron escape time 
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Frequency (Hz) 



Figure 3. Spectrum of the Crab Nebula fitted by our model. The data points are obtained from Baldwin (1971) and Maci'as-Perez 
ct al. (2010) for the radio band; Ney & Stein (1968), Grasdalen (1979), Geen et al. (2004) and Temim et al. (2006) for the infrared; 
Vcron-Cetty & Woltjer (1993) for the optical; Hennessy (1992) for the ultraviolet, Kuiper (2001) for the X-rays and soft 7-rays; and 
Abdo et al. (2010), Aharonian et al. (2004), Aharonian et al. (2006), and Albert et al. (2008) and for the gamma-rays. 



is considered (thus ignoring all other timescales) and Bohm 
diffusion is used. 

Table [2] shows the parameters for each of the models 
needed to obtain a good fit of the Crab Nebula data at the 
current age. The column labelled value are the parameters 
of the complete model of §2 (the origin of each parameter is 
commented in Table [TJ. The dots appear when no change is 
needed from those values. 

Given that the observational parameters such as the 
age, the breaking index, the period, and the period deriva- 
tive are fixed, they continue to determine to and r c in all 
models. For the TDE model, the break energy and the shock 
radius fraction (and, in consequence, the maximum energy 
at injection) have decreased. A smaller shock radius dimin- 
ishes the number of VHE electrons, which is necessary due 
to the lack of energy losses affecting the population, and 
the smaller energy break corrects the lack of radio flux. The 
initial spin-down luminosity is smaller because the lack of 
losses makes electrons' disappearance slower. The magnetic 
field fraction is larger to power the SSC contribution, and 
the energy break increases to compensate the lack of escap- 
ing particles at low energies and correct the radio flux. In 
the ADE-T case, we take the expansion velocity given by 
Tanaka & Takahara (2010) of 1800 km s~\ which gives a 
radius for the PWN of 1.7 pc. This means that the syn- 
chrotron radiation is confined in a smaller volume, so the 
synchrotron photon density is larger and the magnetic en- 
ergy fraction needed to obtain the correct SSC contribution 
is smaller than in our case. Note that the minimum and max- 



Table 2. Comparison of the values used or obtained in the differ- 
ent fits of the Crab Nebula today. Meaning and units of variables 
are as in Table 1. We use dots for those parameters which have 
the same values as in the complete model. 



Symbol 


Value 


ADE-T 


TDE-Z 


L(t age ) 


4.5 X 10 38 




2.5 x 10 38 


7mi„(*) 


1 


10 2 




7t»«.(t) 


7.9 X 10 9 


7 x 10 9 (fixed) 


6.5 x 10° 


7b 


7 x 10 5 


7 x 10 s 


9 x 10 5 


01 


1.5 






«2 


2.5 






e 


1/3 






La 


3.1 X 10 39 




1.7 x 10 39 


B(t aBe ) 


97 




93 


V 


0.012 


0.006 


0.015 


Rpwn 


2.1 


1.7 


1.9 


TCMB 


2.73 






W FIR 


0.25 






T FIR 


70 







WFIR 


0.5 
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imum energy at injection are also fixed in time, according 
with the values used in Tanaka and Takahara (2010). 

It is clear that at the current age, and particularly due 
to the fact of the strong SSC domination of the Crab neb- 
ula, one can find acceptable sets of parameters in both ap- 
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proximated models that fit the data well. However, this 
does not mean that the time evolution of these models 
would be similarly close to the complete analysis. Figure [5] 
and [6] compare the evolution of the results of the complete 
model and the approximate ones, for the electron popula- 
tion and photon spectra, respectively. Differences increase 
with the time elapsed off the normalization age (the cur- 
rent Crab nebula), and are clear at a few hundred and 
a few thousand years. To have a better idea on how the 
spectra are changing, we use our model as reference data 
and calculate the relative theoretical distance of the ADE- 
T and TDE-Z models with respect to it as a function of 
frequency, both for the electron population and spectrum. 
We thus compute the theoretical distance as Distance = 
\complete — approximate\/ complete, so that distance times 
100% is the percentile value of the deviation. These results 
are given in Figure and [8] The dips visible in these figures 
correspond to the crossing of the curves (lepton Lorentz fac- 
tors or photon energies where both the approximate and 
complete models, as plotted in Figures [5] and [6] coincide). 
The recovering of the curves after these dips correspond to 
the use of the absolute value in the Distance definition. Vi- 
sually removing these dips gives an idea of the average devi- 
ation between the approximate and complete models across 
all energies. Note that the peaks in the relative distance evo- 
lution are broader, and represent bona-fide increases in the 
deviation of the approximate results. 

Regarding the underlying electron population, we see 
that differences among models range from 10% to 100% and 
beyond. The TDE-Z model have a more deviating behavior 
than theADE-T models at later ages. Regarding the photon 
spectrum deviation, we find that for ages close to t age = 940 
yr and lower, the relative distance of all models with respect 
to the complete one is below 40% with the exception of the 
frequency range between 10 22 and 10 23 Hz, where there is a 
transition between the synchrotron and IC dominated radi- 
ation. At larger times, deviations can be larger than 100%. 
From 2 to 10 kyr, the relative distance in the optical range 
and in gamma-rays increases with age. As soon as the Crab 
nebula is let to evolve beyond a few thousand years, and 
consistently with the results found for the electron popula- 
tion, the relative distance between the spectra of the com- 
plete and the approximate models goes up to a factor of 
a few (i.e., percentile distance is a factor of a few 100%) 
over large portions of the electromagnetic spectrum. These 
changes in the nebula evolution are only the by-product of 
the approximations used in the models and do not represent 
the expected behavior of the source. 



energy break at ~ 0.35 TeV, a high (low) energy index of 2.5 
(1.5), and a magnetic field of 97 ^iG (in complete agreement 
with morphological studies) fits the data perfectly. Filamen- 
tary structures or flares have not been treated in our model. 
Other losses, associated with neutrino emission or hadronic 
interactions were not taken into account in any of the cases 
presented. 

With the complete model at hand, we have analyzed 
which are the consequences of approximations when models 
ignore losses, photon backgrounds, or escape processes. In 
particular, we analyzed the impact of different approxima- 
tions made at the level of the diffusion-loss equation that 
allows converting it in advective or in a time-derivative-only 
one. Swinging on the parameters one can achieve a relatively 
good fit to the data of the Crab nebula today. However, 
the time-evolution of the electron population and the pho- 
ton spectrum deviation from the complete analysis is larger 
than 100% for these same models, when they evolve in time 
off the normalization age. This puts in evidence the risks of 
considering approximations when studying time evolution, 
as well as, equivalently, when members of a population ob- 
served at different ages are analyzed with the intention of 
extracting statistical conclusions. 
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APPENDIX 

Energy losses 

Synchrotron losses 

In terms of the Lorentz factor of the lepton, we consider the 
synchrotron losses described by 

Wi(7,*) = -| — U B (th 2 , (12) 
3 m e c 

where ar = (87r/3)rQ is the Thomson cross section, ro is 
the electron classical radius, and Us(t) = B 2 (t) /8n is the 
energy density of the magnetic field. We have also assumed 
that the particles are relativistic, i.e., /3 ~ 1. 



5 CONCLUDING REMARKS 

In this work we have introduced a leptonic, time-dependent 
model of PWNe. We have considered the complete time- 
energy dependent diffusion-loss equation to compute the 
lepton population. Full Klein-Nishina cross-section for 
multiple-photon-field inverse Compton, bremsstrahlung, 
synchrotron, and self-synchrotron Compton spectra were 
computed, and their corresponding losses were considered. 
The model has allowed, based on fitting against the injec- 
tion parameters and magnetic field, to reproduce the current 
data for the Crab nebula from radio to TeV. We find that an 



Inverse Compton losses 

We consider the inverse Compton losses described by means 
of the exact expression of the Klein-Nishina cross section, 

3^ i r , 

7/c(7) = -4^^y vfAvf 

x J™ ^/(«, r e )*(l - q)0 (q - -L) dvi, (13) 

where h is the Planck constant, are the initial and final 
frequencies of the scattered photons, 9 is the Heaviside step 
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Figure 5. Electron distribution of the Crab nebula computed for different ages using the complete model, together with the obtained 
results under the ADE-T, and TDE-Z approximations. 



function, and n is the photon background distribution. The 
other terms are defined as usual, 



f(q,T e ) = 2q]nq+ (1 + 2g)(l - q) + gfqfj^, ( 14 ) 

with r e = Arfhvi /m e c 2 , and q = hvf /(T E { r ym e c 2 — hvf)). 
The magnitude F E indicates the regime of the IC energy 
losses. If F e <C 1, the scatter happens in the Thomson limit 
and for r e ^> 1, it happens in the extreme Klein-Nishina 
limit. 



Bremsstrahlung energy losses 

The general expression for the bremsstrahlung losses is 

jBrems = ~NV f k^-dk, (15) 

Jo dk 

where A*' is the number density of particles in the medium, 
v is the velocity of the electrons, k = hv/mc 2 is the pho- 
ton energy in units of the electron rest energy and da/dk 
is the bremsstrahlung differential cross section. The veloc- 
ity v can be expressed in terms of the Lorentz factor as 
v — Cy/y 2 — I/7. We consider two contributions in the 
bremsstrahlung losses: the electron-atom bremsstrahlung 
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Figure 6. Photon spectrum of the Crab nebula computed for different ages using the complete model, together with the obtained results 
under the ADE-T, and TDE-Z approximations. 



due to the interaction of the electron with the electromag- 
netic field produced by the ionized nuclei of the interstel- 
lar medium (ISM) and the electron-electron bremsstrahlung 
produced by the electrons also present in the ISM. This sec- 
ond contribution is the most important and increases with 
energy, but we include the electron-atom bremsstrahlung for 
completeness at lower energies. 



energy range is (Haug 2004) 



f 



7 



— aarZ 2 - 
it 7 + p 



lm( 7 + p)-- 
P 3 



2 2 



19 2 
675 7P 



0.06 



P 



(16) 



In the case of the electron-atom bremsstrahlung, an ac- 
curate approximation for the integral above in the whole 



where a ~ 1/137 is the fine-structure constant, Z is the 
atomic number, and p = \/7 2 — 1 is the linear momentum 
of the electron. The electron-atom bremsstrahlung energy 
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Figure 7. Relative distance of the results for the electron distribution between the complete model, the ADE-T, and TDE-Z approxi- 
mations for the Crab Nebula at different ages. 



losses have the form 
3 



■ e — a 

rems 



7 



aarcS „ 
3 



[7 ln(7 + p) 



f/ 2 2 19 2 
+ e ~ -0.06— 

7° \9 675 7 



with 



5 = ^ ^ 2 iVz = ATh 



(17) 



(18) 



and where Afe is the number density of hydrogen in the 



medium and Nz, the number density of the other elements 
with atomic number Z. 

The electron-electron Bremsstrahlung loss is (Blumen- 
thal and Gould 1970) 



. e-e 3a 



V ZNZ )' 



-(7-1) 



ln(2 7 ) 



(19) 
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Figure 8. Relative distance of the results for the photon spectrum between the complete model, the ADE-T, and TDE-Z approximations 
for the Crab Nebula in different ages. 



Adiabatic energy losses 

The relativistic form of the adiabatic losses is 



lad = 



7- 



(20) 



We consider the PWN as an uniformly expanding sphere, so 
the expansion velocity of the gas can be written as 



v(r) = v PW N(t) 



Rp\VN(t) 



(21) 



Applying the divergence operator in spherical coordinates, 
we get V • v = (l/r 2 )(dv(r)/dr) = 3(v PWN (t)/R PWN (t)), 



and substituting in equation 1201 the adiabatic energy losses 
are 7 a d(7,i) = —(vpwN(t)/RpwN(t))^. We use the free ex- 
panding expansion of the PWN given by van der Swaluw 
(2001). The radius of the PWN is given by 



( -gr - 



V t, 



(22) 



where Vo is the velocity of the front of the ejecta and has 
the form 



Vo 



lOgo 

3M ei ' 



(23) 
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Eq and M e j are the energy of the supernova explosion and 
the mass ejected respectively. The constant C is written as 

-1/5 



c = 



+ 



289 



15(jpwn-1) 240 



(24) 



with 7pwiv = 4/3 since we consider the PWN material as a 
relativistically hot gas. The velocity of expansion can be eas- 
ily obtained doing the derivative of equation (|22p . Applying 
the expressions for RpwN(t) and vpwN(t), we get 



6 7 

T ad--g-, 



(25) 



which differs a factor 6/5 from a ballistic approximation. 

The evolution of the PWN depends also on physical 
parameters of the previous supernova (SN) event, like the 
energy of the explosion and the ejected mass. The age of 
the pulsar has to be less than the Sedov time, which can be 
calculated as (Gaensler & Slane 2006) 



WMr. 



5/6 



E s 



10 51 erg 



-1/2 



"0 



-1/3 



(26) 

where M e j, Esn is the mass and the energy ejected in the 
supernova explosion, and no is the number density of the 
medium. After the Sedov time the PWN is not expanding 
freely due to the interaction with the reverse shock of the 
supernova remnant (SNR) and a dynamic model is needed 
to account for its expansion. 



Photon spectra 

We briefly summary here the expressions that the code uses 
to compute the different luminosities. 



Synchrotron spectrum 

The synchrotron luminosity is (Ginzburg & Syrovatskii 
1965, Blumenthal and Gould 1970) 

f- oc 

L syn {v,t)= N('y,t)P syn (u,'y,B(t))d 1 , (27) 
Jo 

where P SV n{v, 7, B(t)) is the power emitted by one electron 
spiraling in a magnetic field 



P syn (y,~(,B(t)) = 



fc(7,B(*)) 



(28) 



where v c is the critical frequency, F{x) = x J°° K 5 / 3 (y)dy, 
and K 5 / 3 (y) is the modified Bessel function of order 5/3. 
N(-y, i) is calculated solving the diffusion equation explained 
in Section 12.11 The dependence in time points the need to 
recall that the magnetic field has to be computed at the 
same age as the luminosity. 



Inverse Compton spectrum 

The scattered photon spectrum per electron is Blumenthal 
and Gould 1970) 



D , ,n 3a T chu 
P IC { 1 ,v,t) = ~—^ ^ 



n{vi) 



/(?,r e )di/i, 



(29) 



where i>i is the initial frequency of the scattered photon, v is 
the final frequency of the photon after scattering, and n(i>i) 



is the photon target field distribution. The expressions for 
f(q,F e ), F e and q were given above. To obtain the luminos- 
ity, we multiply by the electron distribution and integrate 
in the whole energy range, obtaining 



Lic{v, t) = -archis 



N(rr,t) 



n{i>i) 



/(9,r e )d«/i. 



(30) 



Synchrotron Self-Compton spectrum 



We make use of Eq. (|30[) and the target photon field density 
given by (see Atoyan & Aharonian 1996) 



nssc{v, R ayn (t),t) 



L syn {y,t) U 
4irR1 yn (t)c hv 



(31) 



where R ayn (t) is the radius of the volume where the syn- 
chrotron radiation is produced, and U ~ 2.24 is the mean 
in a spherical volume of the function U(x), given by 



U = 



r «f 
Jo 



RsynW 

— TFT 



x 2 U(x)dx 



f R PW. 

Jo 



777T 



= 3 



r(t) 



Rsyn (^) 



x 2 dx 

Rsyn(t) 

PWN X 2 JJ^ X ^ X ^ 



where 



U[X)= 2J0 



V i x + y 
- In i rdy. 

x \x — y\ 



(32) 



(33) 



The function U(x) was given by Atoyan & Nahapetian 
(1989) to compute the number density of photons at a given 
distance (in this case, R ayn {t)) assuming isotropic emissivity 
of the synchrotron radiation in a spherical source. The value 
of this function at x = is 3 and decreases until 1.5 when 
x = 1. If we consider a photon distribution with a radius 
greater than the radius of the PWN, then U(x) ~ l/x 2 and 
n(u, Rpwn ,t) — L(v,t)/(4TvRp WN c). In this case, the pho- 
ton field would depend on the PWN radius. For simplicity, 
we assume Rsyn/ Rpwn = 1 for all cases. 

Bremsstrahlung spectrum 

The bremsstrahlung luminosity is computed as 



LBrems(v,t) = —aarhcS 

Z7T 



° jV(7j) 



where S is given in Eq. (|18[) . jij are the Lorentz factors 
of the initial and final electron passing through a medium 
containing different kinds of particles and producing photons 
with energy hv. The kinematic condition 7; — 7/ = hu/nieC 2 
fix the final energy of the electron in the integral given an 
initial energy 7; and the energy of the photon produced hv. 
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